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Abstract 

We introduce an algorithm based on a method of snapshots for computing approximate balanced 
truncations for discrete-time, stable, hnear time-periodic systems. By construction, this algorithm is 
applicable to very high-dimensional systems, even with very high-dimensional outputs (or, alternatively, 
very high-dimensional inputs). An example is shown to validate the method. 

I. Introduction 

We consider the model reduction problem for stable, discrete-time periodic systems, using an 
approximation of balanced trunctation. In particular, we combine the lifting approach developed 
in [4] with a snapshot-based approximation described in [15], which is tractable even for very 
high-dimensional systems, on the order of millions of states. 

Several different algorithms are available for extending balanced truncation from linear time- 
invariant systems to time-periodic systems (see, e.g., [13], [9], [4]). Here, our interest is in 
systems with very large state dimension, on the order of tens of thousands or milUons, for which 
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the algorithms used in these previous approaches become intractable. Such systems arise in the 
control of systems governed by partial differential equations, for instance in fluid mechanics, 
where often one approximates the full infinite-dimensional dynamics by a high-dimensional 
discretization. Periodic systems may arise in these settings either as linearizations around a 
periodic orbit (e.g., vortex shedding [14]), or as linearizations of a system with periodic open-loop 
forcing (e.g., periodic pulsed blowing). The goal of this paper is to describe a model reduction 
procedure that closely approximates balanced truncation, yet is computationally tractable even 
for these high-dimensional systems. The resulting reduced-order models may be used for control 
synthesis, or other studies where the full high-dimensional system is unwieldy. 

We suppose that while the number of states and outputs are both very large, the number of 
control inputs is small (by duality, we may alternatively assume that the number of outputs 
is small). This case is also typical in practice, in which one often has a small number of 
actuators, or sources of disturbances, that one wants to model. We also assume that the system 
is asymptotically stable. Key ideas in the method presented here involve the computation of 
the balancing transformation directly from snapshots, without computing the controllability and 
observability Gramians, which become prohibitively large for these systems; and two different 
output projection methods for tractable computation with systems with high-dimensional output 
spaces. 

The paper is organized as follows: in Section|n]we summarize an existing approach to balanced 
truncation for time-periodic systems, and the corresponding lifted time-invariant reformulation. 



In Section III we present the main results, an approximate method based on snapshots, and an 
output projection based on the lifting approach that makes the snapshot method feasible for large 
numbers of outputs. In Section llVl we demonstrate the method on a numerical example. 



II. Background on balanced truncation and linear periodic systems 
A. Balanced truncations for linear time -invariant systems 
Consider a discrete linear time-periodic system 

x{k + 1) = A{k)x{k) + B{k)u{k) 

(1) 

y{k) = C{k)x{k), 
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with state variable x{k) G C", input u{k) G and output y G C^, and k E Z, where the 
matrices A{k) G C"^'^, B{k) G C'^'p and C{k) G C"^'^ are T-periodic (e.g., A{k + T) = A{k) 
for all k). 

We assume that this system is asymptotically stable, in the following sense: defining F{j, i) — 
Aj-iAj^2 ■ ■ ■ Ai for j > i, with F{i, i) — Inxn, the n x n identity matrix, we consider the case 
where the spectral radius p {F{j + T, j)) is smaller than 1. Note that periodicity implies that the 
non-zero eigenvalues of F{j + T,j) are independent of time j [6], whereby asymptotic stability 
is equivalent to uniform geometric decay, as in the LTI case. 

Many different approaches are available for model reduction, including singular perturba- 
tion [8], Hankel norm reduction, and balanced truncation. Balanced truncation [13] has become 
a widely used method for linear systems, since it has a priori error bounds comparable to 
other methods, and is computationally tractable, at least for systems of moderate dimension, say 
n < 10^ However, for very large systems, exact balanced truncation becomes computationally 
intractable, as the procedure involves finding coordinate transformations that simultaneously diag- 
onalize non-sparse nxn matrices (controllability and observability Gramians). To overcome this 
difficulty, a snapshot-based method has been proposed for approximate balanced truncation [15], 
sometimes referred to as balanced proper orthogonal decomposition (balanced POD). Here, the 
idea is to use "empirical Gramians" defined using snapshots of a simulation of the system, as 
in [10], and to compute the balancing transformation directly from the snapshots, without ever 
computing the Gramians. 

The method of balanced truncation has been applied to linear time-periodic systems using 
several different approaches. One approach, used in [9], [11], [18], [16], is to perform T separate 
balanced truncations along a period. An alternative, suggested by [4], is a lifting approach, in 
which balanced truncation is applied to a time-invariant reformulation. However, both of these 
approaches require the solution of Lyapunov equations or inequalities and are not computationally 
tractable for very large systems. Thus, the objective of our work is to obtain more computationally 
efficient algorithms for computing approximate balanced truncations for these large systems. 

Before reviewing the theory for periodic systems, first recall the main idea of balanced 
truncation for Hnear time-invariant (LTI) stable systems (for more detail, see standard textbooks, 
such as [19], [3]). One begins by defining controUabiHty and observability Gramians Wc and Wo 
to measure to what degree each state is excited by an input, and each state excites future outputs. 
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respectively. One then seeks a coordinate transformation $, called balancing transformation, so 
that the transformed Gramians Wc ^ ^~^Wc{^^^)* and Wo ^ ^*Wo^ are equal and diagonal 
(hence balanced): $" Wc($"^)* = $*PVo$ = S = diag(cri, ■ ■ ■ , where ai ^ ■ ■ ■ ^ o-„ ^ 
are the Hankel singular values. Finally, in these new coordinates, one truncates the states that 
are least controllable and observable, corresponding to the smallest Hankel singular values. 

B. Controllability and Observability Gramians for linear periodic systems 

For time-periodic systems, one may also define controllability and observability Gramians, 
as in [18]. The controllability Gramian at time j, Wc{j), and observability Gramian at time j, 
Wo{j), for the complex stable periodic system ([T]) are defined similarly by the following two 
positive-semidefinite n x n matrices 

W,{j):= J2 F{j,t + l)B{t)B{trF{j,z + ir 

(2) 

oo 

WoU) ■.= J2F{z,jrC(zrC(z)F{z,j), 
i=j 

where the superscript * denotes the adjoint of a linear operator. See the Appendix for the specific 
properties these Gramians satisfy. Both Wc{j) and Wo{j) are T-periodic. Also, for each j, Wc{j) 
and Wo{j) satisfy the following respective discrete Lyapunov equations: 

A{j)W,{j)A{jy - W,{j + 1) + Bij)B{jy = 

(3) 

AuyWoU + i)AU) - WoU) + cuycu) = o. 

As discussed in the Appendix, when convergent, the definitions (|2]) (equivalently, ([3])) extend to 
general linear time-varying systems of the form ([T]), with the same interpretations in terms of 
input-to- state and state-to-output mappings. 

C. Traditional approach of balanced truncation for periodic systems 

If the dimension n of the system is not too large, it is possible to solve the above pairs of 
Lyapunov equations for Gramians Wc{j), Wo{j) for each j, 1 ^ j ^ T. One can then do balanced 
truncation at each time j. Refer to [9], [11], [16] for detailed discussions. Note that Farhood 
et al. [4] also compute the Gramians/generalized Gramians by solving corresponding Lyapunov 
equations/inequalities at each time step; however, the balanced truncation is then realized in a 
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lifted time-invariant setting, not in the original periodic one. In this paper we will also use the 
Ufting approach, in conjunction with a method of snapshots as in [15]. 



_y{j + {t + l)T-l)_ 

u{j+tT) 
u{j+tT + l) 

u{j + {t + l)T-l) 



(4) 



D. The lifted system 

As discussed in [5], [4], the linear periodic time- variant system ^ can be equivalently 
rewritten as a linear time invariant (LTI) system: 

Xj{t + 1) = AjXj{t) + BjUj{t) 

ijjit) = CjXj{t) + DjUj{t) 

where j is fixed and 

x,{t) = x{j+tT) 

yU + tT) 

yU + tT + l) 



Febi-uary 5, 2008 



DRAFT 



6 



and 

= F{j + T,j) = Aij +T- l)A{j + T - 2) ■ ■ ■ A{j) 

B,= [E(j + r,j + i)s(j),-- - , 

E{j +T,j+T- + T-2), B{j +T-1)] 

C{j+T-1)F{j+T-1,jI 

o" 

Fj,2,i 

Fj,T,l Fj,T,2 • ■ • 0_ 

where = C(j + i — + i — 1, j ' + k)B{j + A; — 1). We refer to the time invariant 

system (j4]) as a lifted system at time j, in which Aj, Bj, Cj and Dj are constant matrices. The 
lifted system may be viewed as a Poincare map of the original periodic system in the state space 
at times (j + tT), while the input and output information for an entire period is kept at each 
iteration of the Poincare map. 

The controllability Gramian Wjc and observability Gramian Wjo for this LTI lifted system are 
conventionally defined by the two positive-semidefinite n x n matrices 

oo 

w,.:=J2a^b,b; (ij)* 

r (5) 

1=0 

They respectively satisfy the following discrete Lyapunov equations: 

A,W,J* - W,, + B,B* = 

(6) 

A*W,,A,-W,, + C*C,=0. 
The above setting allows one to study the lifted LTI system instead of the original periodic 
one such that the well-developed balanced truncation techniques for LTI systems are available. 
We emphasize that in the lifted LTI system, only one time of balanced truncation is needed. 
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However, the drawbacks of working on the lifted system are that it has even higher dimensional 
inputs and outputs, and that we usually do not have the explicit form of the lifted system in 
computations. The approach taken in this paper is as follows: by using the equivalence between 
the LTI (lifted) and time-periodic systems, we apply the balanced POD procedure in the LTI 
setting but do computations in the original periodic setting. To realize this plan, first we list 
the following statement establishing the equivalent relations between the Gramians of the lifted 
system and the original periodic system. Proofs are by direct calculations, using periodicity. 

Proposition 2.1: Consider an arbitrary time j. The controllability and observability Gramians 
at time j of the periodic system ([T]), i.e., Wc{j) and Wo{j) given by ^ are respectively equal 
to the controllability and observability Gramians Wjc and Wjo of the lifted system ^ at time j 
given by (|5]). 

An immediate corollary from the above statement is that, as stated in [4], [11], solutions 
of the time-varying Lyapunov equations ([3]) are equivalent to the solution of the time-invariant 
Lyapunov equations (|6]) in the lifted setting. This result plays an important role in traditional 
hfting approaches of balanced truncations for periodic systems, where the Gramians are found 



by solving Lyapunov equations. In the following section, we show how Proposition 2.1 allows 
us to apply a snapshot-based algorithm for computing balancing transformations for periodic 
systems, thereby avoiding the need to solve Lyapunov equations. 

in. Balanced POD for linear periodic systems 
A. Method of snapshots for computation of empirical Gramians for time-varying discrete systems 

The method of snapshots provides a numerical approximation of Gramians directly based on 
their definitions. Refer to [15] and the references therein for the application of this method to 
hnear time-invariant systems. 

For time-varying discrete systems, first, we define the empirical controllability Gramian and 

empirical observability Gramian at time j for the periodic system ([T]) as 

i-i 

W,e{3\m):= J2 F{j,z + l)B{i)B{i)*F{j,z + iy (7) 

i=j—Tn 
j+m—l 

Woe{r,m):= J2 F{i,jrC{i)*C{i)F{i,j), (8) 
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where m is a nonnegative integer. Clearly, limm_^oo W^ce(j; = Wc{j); limm~*oo Woe{j]rn) = 
Wo{j). Further, the following result gives upper error bounds of ||M/c(j) — W^ce(j; ^^i) || and 
\\Wo{j)-Woe{j;m)\\. 

Lemma 3.1: Consider an arbitrary time j. Let m = IT, where / is a nonnegative integer. For 
the linear periodic system ([T]), 



mm 

\\Wo{j)-Woij;m) 



^ \\F{j + T,jyr 
^ \\F{j + T,jyf, 



(9) 



wwoim 

where an induced norm is used. 

The proof is based on the following result for linear time-invariant systems. 
Lemma 3.2: Consider a linear stable time-invariant system 

x{k + 1) = Ax{k) + Bu{k) 
y{k) = Cx{k), 

where A, B and C are constant matrices and whose controllability and observability Gramians 
are Wc = YZo^'BB* {A')* and Wo = E^o (^0* C*CA\ Define the empirical controllability 
and observability Gramians for this LTI system as 



i=0 

l-l 



Woe{l) := J2 {^TC*CA' 



(10) 
(11) 



i=0 



where / is a nonnegative integer. Thus, 



\Wr.-Wr,. 



\Wr. 



€ \\A 



l\\2 



\Wo\ 



€ \\A 



l\\2 



where an induced norm is used. 

Proof: Under any induced norm, 

WWr-Wr., 



J^A'BB* [A'Y 



i=i 



^ P llllW^cll 
Similarly, we have \\Wo - Woe{l)\\ ^ WA^fWWol 



A'WciA'y 



lA'fWWrW. 
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Proof: [Proof of Lemma 3.1 1 By Proposition 2.1 Wc{j) = Wjc- Also, it is straightforward 



to show Wceij',m) = Wjce{l), where Wjce{l) is the empirical controllability Gramian of the 



lifted LTI system at time j, defined in form of (10). Thus 



and the result follows by applying Lemma 3.2 to the lifted LTI system at time j. The result for 
the observability Gramian follows from a similar argument. ■ 
Indeed, for any e > 0, there is an induced norm ||-|| such that ||F((j+T, j) || ^ p ' + T, j))+e 
(by Lemma 5.6.10 in [7]). Thus, with an e satisfying p {F{j + T, j)) +e < 1, Lemma XT] implies 
that under a certain matrix norm, 

l|waj)-M'=a;jnU<,,,f,, + r.i)) + .f 



where m = IT. Thus, a large enough m guarantees that the error between empirical and exact 
Gramians is small. This result encourages one to use empirical Gramians instead of exact ones 
to realize approximate balanced truncations, whenever it is difficult to compute exact Gramians. 
The following subsections introduce the method of snapshots by which we can calculate those 
empirical Gramians for linear periodic systems. Also note that it is clear that the condition 
m = ZT is not crucial to obtain a result similar to that in the Lemma 13.11 

1) Computation of the empirical controllability Gramian: Without loss of generality, let 1 ^ 
j ^ T. The empirical controllability Gramian Wce{j]m) can be obtained by running a series of 
impulse-response simulations for system ([T]) as follows. 

Consider the d-th column of B{k), denoted by [B{k)Y. To obtain a finite sum J2i=j-m -^0' ^ + 
l)[B{k)]''- (^[B{k)]'^y F{j,i + 1)*, one starts the first simulation with the initial condition x{j — 
m) = and a unit impulse u{j — m) = [0, ■ ■ ■ , 0, 1, 0, ■ ■ ■ ,0]^ whose d-th entry is 1. By running 
m steps, we obtain 

x{j -m + 1) = [B{j - m)f; 

x{j -m + 2) = F{] - m + 2,j-m + l)[B{j - m)Y; 
xU) = F{j,j -m + l)[B(j - m)Y. 
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Notice that x{j) is a snapshot we need since x{j)x{j)* appears in the finite sum. Similarly, we 
run the second simulation with initial condition x{j — m+l) = and unit impulse u{j — m+l) — 
[0, • • • , 0, 1, 0, • • • , 0]^, and by running (m — 1) steps we get x{j — m + 2) = [B[j — m + 
■ ■ ■ , x{j) — F{j,j — m + '2)[B{j — m + 1)]*, where the x{j) will be used to construct the finite 
sum. Repeat the above process until the m-th simulation, with initial condition x{j — 1) = and 
unit impulse u{j — 1) = [0, ■ ■ , 0, 1, 0, ■ ■ ■ ,0]^, where only one step of simulation is needed to 
obtain x{j) = [B{j — 1)^. By running the m simulations with m + ■ — hi = m(m + l)/2 steps 
in all, we can thus construct the sum I]-=i-m ^(i' ^ + ^)[B{k)Y {[Bik)^)* F{j,i + 1)*. Since 
B{k) has p columns, we need to run mp simulations, with m{m + l)p/2 steps total, to obtain 
mp snapshots. Clearly, we have 

Proposition 3.3: (Controllability Gramian from snapshots) Define an n x mp dimensional 
matrix 

X(j;m)^[F{j,j-m + l)[B(j-m)]\--- ,[B{j-l)]\.-. , 

(12) 

F{j,j -m + l)[B{j - m)Y, ■■■ ,[B{j- l)]"] , 
whose columns are the snapshots obtained through the mp impulse-response simulations de- 
scribed above. Then 

X{j;m)X{j;m)*^WUj;m). (13) 

Note that physically the empirical controllability Gramian is based on considering the total 
influence of the past history of impulse inputs on the 'current state' x{j) by neglecting those 
impulses u{k), k < j — m, whose influence on the current state has mostly decayed with 
sufficiently large m. 

The above procedure is generally valid for calculation of the empirical controllability Gramian 
of a time-varying stable system. The T-periodic feature of our system can substantially save the 
computational effort. For instance, the snapshot F{j, j — m + T + l)[B{j — m + T)^ is obtained 
by running the '(T + l)-th' simulation mentioned above for m — T steps, with initial condition 
x{j — m + T) — and corresponding unit impulse. However, to run this simulation is not 
necessary, since F(j,j -m + T + l)[B{j -m + T)]"^ = F(j - T,j -m + l)[B{j - m)]'^, 
which has been calculated in the 'first' simulation we mentioned above at the (m — T)-th step. 
Indeed, it is clear that for a T-periodic system, when m > T, to obtain the mp snapshots 
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for construction of X{j;m), one needs only run Tp impulse-response simulations: For each 
column [B{k)Y, starting from x{j — m) — run a simulation for m steps. Then, start from 
x{j — m + 1) — 0, running for m — 1 steps. Repeat the process until the simulation starting from 
x{j — m + T — l) — and running for m — T + 1 steps. The data set obtained is sufficient for 
construction of X{j;m). 

2) Computation of the empirical observability Gramian: To calculate Woeij'jin) for a fixed 
j, 1 ^ j ^ T, by the method of snapshots, one needs to construct an adjoint system 

z{k + l)^A{k)z{k) + C{k)v{k) (14) 

where k — j, - ■ ■ , j + m — 1, the state z{k) e C*, the control input v{k) e C^, and 

A{k) = A{2j + m-k-l)\ C{k) = C{2j + m-k- 1)*. 

As above, here we run a series of impulse-response simulations for the adjoint system. Consider 
the d-th column of C(k), [C(k)]'^ = [C{2j + m-k- One starts the first simulation with 
the initial condition z{j) —0 and a unit impulse v{j) — [0, • • • , 0, 1, 0, • • • ,0]^ whose d-th entry 
is 1. By running the simulation for m steps, we obtain 

z{j + l)^[C{j + m-iy]'; 

z{j + 2)=F{j+m-l,j+m- 2)* [C{j +m- lyf ; 

z{j +m)= F{j +m- 1, jT [C{j +m- 1)*]" , 

where z{j + m) is a snapshot we need for computing the Woe{j',rn). The second simulation 
starts with initial condition z{j + 1) = and unit impulse v{j + 1) = [0, • • • , 0, 1, 0, • • • , 0]^, 
and by running (m — 1) steps we stop at z{j + m) = F{j + m — 2,j)* [C {j + m — 2)*]'^, another 
snapshot we need. Repeat the above process until the m-th simulation, with initial condition 
z{j + m — 1) = and unit impulse v{j + m — 1) = [0, ■ ■ ■ , 0, 1, 0, ■ ■ ■ , 0]^, where only one step 
of simulation is needed to obtain z{j + m) — [C{j)*]'^. Since J{k) has q columns, we need to 
run mq simulations, with m{m + l)q/2 steps total, to obtain mq snapshots, with which we have 
the following statement: 
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Proposition 3.4: (Observability Gramian from snapshots) Define an n x mq dimensional 
matrix 

Y{r,m)= \F{j+m-l,jr[C{j+m-ir]\--- ,[C{jr]\--- , 

^ , (15) 

Fij + m - 1, jT [CU + m - l)*]"^ , ■ ■ ■ , [CUYY \ , 

whose columns are the nioq snapshots obtained through the mq impulse-response simulations 
described above. Then 

Y{j;m)Y{j;mr = Woe{j;m). (16) 



Again, the above procedure is valid for general stable time-varying systems. If the system is 
T-periodic, then as for the empirical controllability Gramian case, one only needs to run Tq 
simulations to obtain the mq snapshots. 

B. Balanced truncation using the method of snapshots 

Suppose we have obtained the factors X{j; mc) and mo) mentioned above, where mc, mo G 



Its 



N may be different. Then for the lifted LTI system at time j, by Proposition 2.1 , 3.3 and 3.4 
Gramians are approximated by 

Wjo ^ Woeij] mo) = X{]- mo)X{j; mo)* 

(17) 

Wjo ^ Woe{j; mo) = Y{j; mo)Y{j; mo)*. 
For LTI systems, the method of snapshots presented in [15] gives an algorithm for computing 
the transformation that exactly balances the empirical Gramians Wee and Woe, directly from the 
factors X{j]me) and Y(j;mo), without computing the Gramians themselves: 

Theorem 3.5: (Balanced truncation using the method of snapshots [15]) Let S G C"^" be 
a real diagonal matrix including the non-zero Hankel singular values, obtained by singular value 
decomposition (SVD) of the matrix Y (j ; mo)* X {j ; me) 

Y{j;mo)*X{j;me) = UJ:V*, (18) 

in which a is the rank of Y{j;mo)*X{j;me), and U G C"""^", V G C'"^^^'' satisfy U*U = 
V*V = I ax a- The balancing transformation is then found by computing matrices 

$ = X{j; m,) ^ ^ Y{j; mo)UT.-^'\ (19) 
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If a = n, then $ is the balancing transformation and is its inverse. If a < n, then the columns 
of $ form the first a columns of the balancing transformation, and the rows of ^* form the first 
a rows of the inverse transformation. 

We emphasize that computing the balancing transformation by the method of snapshots 
exactly balances the empirical Gramians. The only approximation here is to use X and Y to 



approximately construct the Gramians as in (17). 



Balanced truncation of order r is then done as follows. Let $i denote the first r columns 
of $, and the first r rows of The reduced state variable Zj(t) E C satisfies Zj{t) = 
'^lxj{t) = \E'ix(j + tT). The reduced model, in the lifted setting, is then given by 

Zj{t + 1) = ^lAj^iz,{t) + ^IBjUjit); (20) 
y,{t) = C,^iz,{t) + DjUj{t), (21) 



In simulations, the reduced output equation (21 1 shall be un-lifted to the original periodic setting: 
For each i, i = 1, ■ ■ - T, 

y{3 +tT + i-l) = C{j+i- l)F{j + I - 1, J>i5i(i) 

T (22) 

+ Y^Dji^i^k)u{j +tT + k-l) 

k=l 

where Dj{i,k) denotes the entry of Dj at i-th row and k-th column. 

C. Output projection method 

When the number of outputs q is very large, direct construction of F(j;mo) as described in 



Section III-A.2 will be computationally intractable, because Tq adjoint simulations are required. 
To overcome this, an output projection method [15] may be used, by which one can substantially 
reduce the number of adjoint simulations. The starting point of this method is to define an 
optimization problem that minimizes the error between the input-output behavior of the original 
system and that of a projected system with a smaller-dimensional output space. For a time- 
periodic system, this optimization problem is cumbersome to define directly. We instead consider 
the lifted LTI system, for which we design a version of output projection, and then relate the 
method back to the periodic system. 

Fix a time j. The lifted system at time j, though with an even higher dimension of output Tq, 
is a standard LTI system, and its input-output behavior can be measured by a sequence of Tq x Tp 
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dimensional impulse-response matrices {Gj{t)}, the i-th column of each Gj{t) representing the 
output response yj{t) corresponding to a unit impulse input ttj(O) = [0, ■ ■ ■ , 0, 1, ■ ■ ■ , 0]^ whose 
i-th entry is 1. By output projection [15], we mean the search of an orthogonal projection on 
C^i with rank f^p < Tq, i.e. a Pj = 6^6* where Qj G C^'^^''°^ 6*6^ = /f^^xfop, that leads to 
a projected lifted system at time j 

Xj{t + 1) = AjXj{t) + BjUj{t)] (23) 

y,{t)p = Pj {C,x,it) + b,u,{t)) , (24) 

which is an approximation of the original lifted system, such that only fop adjoint simulations 
of the corresponding adjoint system are needed for calculation of the empirical observability 
Gramian of the projected system. More details will be discussed soon. The Pj shall be chosen 
such that the input-output behavior of the projected system is as close as possible to that of the 
original LTI system. More precisely, Pj is the solution of the optimization problem 

min I V||G,(t)-P,G,(t)||M (25) 

with respect to some norm on matrices, where Vr denotes a space of rank-r orthogonal projec- 
tions. If we use an induced norm, such as the Frobenius norm 1 1 ■ | |f induced by the inner product 
{A, B) = Tiace(A* B), then the minimization problem has a standard solution: Pj = 0j6j, 
where Qj = 0], ■ ■ ■ , &j°'' in which the column vectors {6* } are the orthonormal eigenvectors 
of 7^ = Xli^o It ^ typical eigenvalue problem in POD reduction and can be 

numerically solved by the method of snapshots [17], where the snapshots are provided by the 
data sets {G'j(i)}*^g obtained through simulations. 

The above output projection generally produces a full matrix Pj. Thus the components of the 
projected yj{t)p = Pjyj{t) no longer cleanly correspond to the outputs of the periodic system at 
different time steps respectively, but the intermixed combinations of them, which is not desirable 
both for physical understanding of the system, and for numerical simulation purposes (recall that 
in simulations we do not really compute with the lifted system, so the projected lifted system 
should be 'unlifted' back to a periodic system for computation). To deal with this difficulty, we 
propose a modified version of output projection, in which we seek a sub-optimal solution that 



solves (25) with a constraint imposing that the orthogonal projection Pj with rank Top takes a 
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diagonal form 



= diag P,(l) 



(26) 



where each q x q block is an orthogonal projection on with rank in form of Pj{i) 
Qj(i)Qj{iy where e C«"'''°f and Qj(iyQj{i) = Iro^xrop- Note that in this case 



e, = diag e,(i),-- - ,e,(T) 



(27) 



Here we need Top = ropT. 



The diagonal Pj makes it easy to unlift the projected lifted system (23) & (24) to a projected 
time-periodic system 

x{k + l) = A{k)x{k) + B{k)u{ky, (28) 
y{k)p = P{k)C{k)x{k) (29) 

for A; = j, ■ ■ • , where the T-periodic orthogonal projection P of rank Top on is defined by 

P{j +tT + z) = P{j + ^) := + 1), 2 = 0, ■ ■ ■ , T - 1, 

(30) 

W7 = 0, I,--- . 

We write P{k) = Q{k)e{k)* where the T-periodic 6 G C^^^'^f is defined by e(j + tT + i) = 
e(j+0 := Q,{^ + l). 

By the above unlifting, though solving for the sub-optimal Pj is no longer a standard POD 
reduction problem in the lifted setting, we can attack it in the periodic setting such that standard 
POD reduction techniques are available. 

To realize that, first we rewrite Gj{t) as 

'G,(t)i 



(31) 



where each block Gj{t)i is a g x Tp matrix. Let ^ 6 ^ (T — 1), 1 ^ c ^ p. The (bp + c)- 
th column of Gj{t)i represents the output response of the original periodic system at time 
{j + tT + i — 1) to the unit impulse u{j + h) = [0, ■ ■ ■ , 0, 1, 0, ■ ■ ■ , 0]^ whose c-th component 
is 1. Thus, for the periodic system ([!]), we define 



G{j+tT + i,j) :=G,(t),+i 



(32) 
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as its impulse-response matrix at time (j + tT + i), since it includes all different responses at 
the current time respectively to corresponding unit impulse inputs during the whole time period 
[j, {j + T — 1)]. This definition matches that proposed in Bamieh and Pearson 1992 [1]. 

The following statement links the constrained optimization problem given above for the lifted 
system to an equivalent optimization problem defined in the periodic system. 



Proposition 3.6: Under the Frobenius norm, the optimization problem (25), with its solution 



Pj constrained in form of (26), is equivalent to the optimization problem 

^T-l oo 



mm 

{P{j +i)&r rop, 

i=o,-,r-i} 



G(j+tT + z,j) 



J=0 t=0 



(33) 



which implies, for each i = 0, 



mm 

{P{j + i)(^rrop} 



-P{j+t)G{j+tT + t,j) 
T — 1, P(j + i) is the solution of the optimization problem 

oo 



.t=o 



(34) 



-P{j+i)G{j+tT + z,j: 



Proof: By direct calculation, using ( |26| ), ( |3T| ), ( |30| ), ( [32] ) and the linearity of the Trace 
operation. ■ 



This statement allows us to obtain the Pj in form of (26) by solving T unconstrained 
optimization problems forP(A;) = Q(k)Q{k)*, k = j, - ■ ■ ,j + T — 1, in the periodic setting, each 
of which is a typical eigenvalue problem in POD reduction with solutions satisfying that the 
Top columns of each Q{k), {0{kyyi^^, are orthonormal eigenvectors of Tl{k) = Ylt^o^i^'^ + 



k,j)G{tT + kjy, i.e., they satisfy 



R{k)e{ky = Xie{ky. 



(35) 



Numerically, the eigenvectors can be solved by the method of snapshots, in which T SVDs 
are needed to obtain 6(j), • • • , 0(j + T — 1). See details of this method in, for example, [17] 
and [15]. The snapshots are the columns of impulse-response matrices {G{j + tT + i, j)}f^Q for 
2 = O,-- - ,T - 1. 

A convenient computational feature of the output projection method is that by periodicity we 
see all the snapshots have already been obtained during the computation of X{j; rric) described 
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in Section III-A.l as long as rric ^ (s + 1)T, except that the data obtained there needs to be 
left-multiplied with a corresponding C matrix. For instance, the matrix C{j)X{j;mc) includes 
all the columns of matrices G{j + tT,j) for w = 1, - ■ ■ , rric/T. 



Suppose we have found the Pj in form of (26). By definition, the observability Gramian of 
the projected lifted system at time j is 



and the observability Gramian at time j for the projected periodic system is 

oo 

WopU) = Y.F{i,jrC{i)*Q{i)Q{i)*C{i)F{i,j). 



(36) 



(37) 



1=] 



The following statement for the projected systems is an analog to Proposition 2.1 for the full 
case. Its proof is by direct calculation. 

Proposition 3.7: The observability Gramian for the projected lifted system at time j, Wjop, 
is equal to the observability Gramian at time j for the projected periodic system, Wop{j). 

Therefore, Wjop = Wop{j) ~ WoPeij'-,'mo), where WoPeij^rrio) is the empirical observability 
Gramian at time j of the projected periodic system obtained by the method of snapshots 
introduced in Section 3.1. Now the input is only rop dimensional, rop <^ q, in the corresponding 
projected adjoint time-periodic system 

z{k + 1) = A{k)z{k) + Cp{k)vr^^{k) (38) 
where k = j, - ■ ■ , j + mo — 1, the state z{k) e C", the control input Vr^p{k) G C'^°p, and 

A{k) = A{2j + mo-k- 1)*, 
Cp{k) = C{2j + mo-k- l)*e(2j + mo-k-l). 

Thus, only Tvop adjoint simulations in total is needed for computing Weop{j] rno). 



Remark 3.1: Instead of seeking sub-optimal solutions for ( [25) ) in form of ( [26| ), ( p7) ), where 
each = 6(j + i — 1), i = 1, ■ ■ ■ , T, are different and given by the solutions of (35), one 

can also alternatively impose a stronger constraint that all the projection matrices at each time 
step are the same, which means 6(j) = 6(j + 1) = ■ ■ • = 0(j + T — 1) := 6, and 



e, = diag[e,--- ,6]. 



(39) 
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It is straightforward to show that, similar to (35), the single 6 is obtained by solving the 
eigenvalue problem 

Re^ = Xie\ (40) 

where {^'j^^i are the columns of 6, and R = EitJ"^ YZoG{tT + k,j)G{tT + k,])*. The 
6 obtained by this setting reflects the overall impulse-responses of the periodic system, mixing 
those at different time steps along each time period. The stronger constraint given here let us 
expect that the solution is even less optimal. However, a numerical advantage is that only one 
SVD is needed in the method of snapshots for solving the single 0. 

D. Summary: procedures of balanced POD 

Following the terminology in [15], the approximate balanced truncation based on the method 
of snapshots given in Section III-A & III-B[ plus the output projection introduced in Section III 
|C[ is named balanced POD. We summarize the computational procedures of the balanced POD 
method for linear asymptotically stable time-periodic systems: 

• Step 0: Pick a time j, 1 ^ j ^ T, as the "base point," based on which the balanced 
truncation will be done. 

• Step 1: Run Tp impulse-response simulations to obtain rricP snapshots and form n x rricP 



dimensional X{j;mc) as described in Section III-A 
Step 2: For each (j + i = 0, ■ ■ ■ ,T — 1 (corresponding to one whole period), left-multiply 
the state responses stored during computing X(j;mc) by corresponding C matrices, and 



then use the method of snapshots to solve the T eigenvalue problems defined by (35) for 
the T output projection matrices P(j + i) along one time period. 



Step 3: Construct the "projected adjoint system" (38) and run Tvop impulse-response sim 



ulations for the adjoint system to form n x moTop dimensional Y{j]mo), as described in 
Section UlTAl 



Step 4: Compute the SVD of ■mo)*X{j; rric) defined in ( 18 ) and compute the balancing 



POD modes for the lifted system given by (19). 



Step 5: Obtain the reduced lifted system in form of (20) & (21 ). For computational purposes. 



rewrite the output equation (21) as (22). 
If the dimension of output q is small, then we can skip Step 2 for output projection, and 



directly run Tq impulse-response simulations for the adjoint system (14) as in Section 3.1 to 
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form the n x nioq dimensional Y{j;mo). On the other hand, in Step 2, instead of solving for 
the time varying T-periodic output projection matrices, one can alternatively seek a single time- 



invariant projection matrix by solving the eigenvalue problem (40) as we mentioned in Remark 
Ol 

Remark 3.2: Note that to pick up an optimal base time j in the sense that the upper error 
bound of the reduced system of a fixed order is minimized, one needs to solve for the Hankel 
singular value matrices S corresponding to the lifted systems at different times during one whole 
period (see Farhood et al. [4]). This method is therefore not attractive or even computationally 
feasible for large systems. 

Remark 3.3: A "dual" version of the above algorithm is readily available for balanced trun- 
cations of linear periodic systems with high-dimensional states and inputs, but only few outputs. 
In that case, one can start with the construction of Y{j]mo) by running a series of impulse- 
response simulations for the adjoint system, since the dimension of outputs is small. Then an 
input projection based on POD reduction shall be done in the same spirit as that underlies the 
output projection to obtain a projected system whose dimension of inputs is reasonably small. 
One then runs a corresponding number of impulse-response simulations for the projected system 
to construct X(j;mc). 

IV. Example 

To validate and demonstrate the balanced POD algorithm, we consider the following example. 
Consider a linear periodic system Q with period T = 5, state dimension n = 30, output 
dimension q = 30, control input dimension p = 1, and {A(k)}l^^ are randomly generated 
diagonal matrices with diagonal entries bounded in [0.16,0.96], guaranteeing the asymptotical 
stability. {B{k)} and {C{k)} matrices are also randomly generated, with entries bounded in 
[0, 1]. The setting is similar to that used in the example in Farhood, et. al [4]. It is not a high- 
dimensional problem, such that exact balanced realizations for the corresponding lifted system 
can be done by solving Lyapunov equations for Gramians of the lifted system, and the result 
can thus be compared with that by the balanced POD approach. 

We pick the "base time" j = I. Choose = rrio = 2T = 10, and test different cases of 
balanced truncation, including balanced PODs with the order of output projection for the periodic 
system at each time step rop = 1, 2, 6 and 10 respectively. Recall that, to the lifted system, the 
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order of output projection is = VopT. Figure [T] shows the Hankel singular values in the 
lifted setting, obtained by the exact balanced realization, the balanced truncation based on the 



method of snapshots given in Section III-B but without output projection, and the balanced POD 
(with T different output projection matrices along one period for the periodic system). Figure 
[2] shows the error plots of the infinity norm, \\G — Gf ||oo/||G||oo versus f, for those different 
cases. Here Gf is the impulse-response matrix of the reduced lifted system of order f. We do 
not show the balanced POD with Top = 10 case since the result is almost identical to that by 
balanced truncation based on the method of snapshots but without output projection. We see 
that the balanced truncation based on the method of snapshots gives a good approximation of 
the exact balanced truncation, and further, the balanced POD, even with low orders of output 
projection rop, generates satisfying results. Figure |3] shows comparisons between balanced POD 
results with the same order of output projection, one set based on T different projection matrices 
along one period in the periodic setting, and the other only having one time-invariant projection 



matrix (see Section III-C). For the cases where Top are low, these two approaches give almost 
identical results, or even the latter one gives better results. However, when the order of output 
projection Top increases, such as r^p = 6, the results based on T different projection matrices at 
each time step are better than those by a single projection matrix, as we expect. 

V. Conclusion and future directions 

We have proposed a version of the balanced POD method to realize approximate balanced trun- 
cation for linear asymptotically stable periodic systems, especially with very high-dimensional 
states and outputs but a small number of inputs. It is a generalization of the balanced POD 
method for linear time invariant system developed in [15]. The development of this balanced 
POD is based on a lifting approach, and key parts include the method of snapshots for computing 
empirical Gramians, and a version of output projection that gives a periodic projected system 
based on POD reduction, with which the number of adjoint simulations needed for computing 
the empirical observability Gramian is substantially decreased. Simulation results given in the 
previous section validate the approach. This snapshot-based approach is also readily applicable 
to high-dimensional systems with high-dimensional inputs and few outputs, where an input 
projection is needed. 

A future direction of this work is to apply the balanced POD method to construct reduced-order 
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Fig. 1. The Hankel singular values aj: exact balanced truncation(n); balanced truncation by the method of snapshots but 
without output projection(o); balanced POD with Vop — 1 (v); balanced POD with rop — 2 (*); balanced POD with Vop — 6 
(+); balanced POD with Top = 10 (x). 

models of high-dimensional (linearized) periodic systems arising in engineering applications and 
then, based on the low-dimensional models, to design closed-loop control laws based on these 
models. For instance, such periodic orbits may arise as periodic shedding in the wake of a bluff 
body [14], or from open-loop forcing at a prescribed frequency, as in the recent results of [12], 
which show that a periodic blowing and suction of flow added at the walls of a channel flow 
may reduce drag. The dimension of states of such systems, including the three components of 
velocity and pressure at each grid point in the channel, can be on the order of 10^, and thus 
reduced-order models of such systems are quite valuable for designing model-based control laws. 

Appendix 

a) Observability Gramian for linear time-varying systems: The observability Gramian 
provides a measure of the influence of an initial state x{j) on future outputs with zero control 
inputs. To see that, first, for a fixed time j, we define a linear operator \l'oj : C" oo) for 

system ([T]) to describe the state-output behavior y = \l'ojx(j) with zero inputs and an initial state 
x{j). More precisely, y{j + k) = C{j + k)F(j + k,j)x{j). A; = 0, 1, ■ ■ ■ . To measure to what 
degree the state x{j) excites the output y, it is natural to compute the square of the induced 
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Fig. 2. Error \\G — Gf\\.x/\\G\\cx,, for exact balanced truncation(n), balanced truncation by the method of snapshots but 
without output projection(o), balanced POD with Top = 1 ( ), balanced POD with Vop — 2 (*), balanced POD with Vop = 6 
(+), and the lower bound for any model reduction scheme (— ). 



norm 

\\y\\l = {y,y)p 

where (■, ■) denotes an inner product, with subscripts specifying the vector space where the inner 
product is defined. The observabiUty Gramian given in (|2]) has the following property: 
Proposition 1.1: 

where : oo) C" is the adjoint of "^oj- And, 

I blip = {xU),WoU)x{j))cn ■ 

Proof: First, it is clear that 

Mil = {y,y)i^ = (^oix(j),^oix(j))p = {x{j),%j^oMj))c« ■ 
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Fig. 3. Time varying T-periodic output projections versus time-invariant output projections: Error j|G — G,-] |oo/| IG] |oo, for 
balanced POD witii Top = 1 ( balanced POD with Top = 2 (*) and balanced POD with Top = 6 (+). Solid lines correspond 
to cases using T different projection matrices along one period for the periodic system, and dashed lines using one single 
projection matrix. The black solid line is the lower bound for any model reduction scheme (— ). 



To explicitly express \E'*j, we consider, for an arbitrary z E oo), 

oo 

= J2{z{z),C{t)F{z,j)x{j))^„ 
i=j 

= /f2F{^,JrC{^rz{^),x{J)\ 

where matrices are considered as linear operators and (■)* denotes the corresponding adjoint. 
So = ET=,FihjrC{zyz{z), and K.i^oMj)) = YZ, F{i, 3yC{iYC{i)F{i, 3)xU) = 
Wo{j)x{j). U 
Note that the observability Gramian has non-negative eigenvalues, the larger ones correspond- 
ing to the more observable states. 

b) Controllability Gramian for linear time-varying systems: Similarly, the controllability 
Gramian provides a measure of the influence of a sequence of input history on the current 
state. In other words, to reach a given current state (if possible), the controllability Gramian 
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measures how much effort from control inputs is needed. For a fixed current time j, define a 
linear operator \E'cj : ^^(— oo, j — 1] ^ C" such that x{j) = "^cjU = Xlt-oo -^0' + ^)B{i)u{i). 
This operator describes the input-state behavior with initial state x(— oo) = and a sequence 
of inputs {u{i)ys^. Consider the 'energy' of input needed for reaching the current state x{j) 
(suppose the system is controllable at time j) 

= {u,u)i2. 

We have: 

Proposition 1.2: The controllability Gramian Wc{j) defined in ^ satisfies 
where : C" ^ P{-oc,j - 1] is the adjoint of \E'cj- And, 

Mil = {xij),w,ij)-Hj))c^ 

where (•)^^ is the inverse of (•). 
Proof: First, 

where one uses the fact that (^;/)* = (^^) . 

Similar to the observability Gramian case, by calculations under standard inner products of 
and /2(-oo,j-l], one obtains (^^2) (i) = B{iy F{],i + l)* z, i = -00, ■ ■ ■ z e C". 

It follows from definition that WdJ) = ^cj^cj- ■ 

Note that the eigenvalues of the controllability Gramian are non-negative and the larger ones 
correspond to the more controllable states. 
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